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Abstract 


The first part of the present thesis describes an experimental technique for the inverse deter- 
mination of thermal conductivity in a 2D square domain based on a finite-difference based 
numerical algorithm, which uses transient temperature data as the input. The algorithm is 
tested for a known conductivity material such as mild steel.lt is found that the algorithm 
is highly sensitive to the measurement errors in the input data, a fact not unexpected in in- 
verse problems which are essentially ill-posed in nature. In addition to the above mentioned 
experimental work, a detailed stability analysis is also performed in order to find the maxi- 
mum permissible measurement error for such problems using simulated temperature profiles 
containing random errors. Interestingly, the study reveals that the maximum permissible 
measurement error is high at the early time and decreases exponentially with time. This 
means that early time temperature profiles should be used as input data. However, even at 
early times maximum allowable error in the temperature measurement is quite small which 
explains why. the present experimental input data did not give correct conductivity value for 
mild steel.Future effort is needed in the direction of improving the experimental set-up so 
that high precision temperature measurement is possible. 
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Greek Symbols 


a relaxation factor 
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Superscripts 

p present time 
p + 1 future time 
' first derivative 

" second derivative, also per unit area 


Special Symbol 
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Chapter 1 

Introduction and Literature Review 


1.1 Introduction 

The use of inverse heat conduction techniques for the determination of thermal prop- 
erties such as thermal conductivity and heat capacity of solids or the estimation of surface 
condition such as temperature and heat flux by utilizing the transient temperature measure- 
ments taken within the medium, has numerous practical applications. 

In the past decade, rapid advancement has been made in the field of materials research. 
This is primarily due to the demand for advanced material with low weight or the ability 
to withstand high temperatures required in the nuclear, electronics and aerospace industries. 
Precise knowledge of thermophysical properties for these materials is essential in many ther- 
mal systems. Specifically, an accurate prediction of thermal conductivity is imperative to 
achieve thermal control system. 

The direct measurement of heat flux at the surface of a wall subjected to a fire, at 
the outer surface of a re-entry vehicle or at the inside surface of a combustion chamber 
is extremely difficult. In such situations, the inverse method of analysis, using transient 
temperature measurements taken within the medium can be applied for the estimation of 
such quantities. 

An inverse problem is much more difficult to solve than the direct one. The reason for 
this is that it is usually ill-posed, i.e., it is very sensitive to measurement errors. An excellent 
discussion of difficulties encountered in inverse analysis is well documented in Beck et al. 
(1985). 

The present work involves the estimation of thermal conductivity by an inverse tech- 
nique using finite-difference. Therefore, the literature review that follows concerns only the 
papers and texts associated with thermal conductivity measurements. 
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1.2 Literature Review 


To date,various methods have been developed for analysis of the inverse heat conduc- 
tion problem involving the estimation of thermal conductivity from measured temperatures 
inside the material. Beck and Arnold(1977) determined the thermal conductivity by mini- 
mizing the errors between the measured and calculated temperatures in a least squares sense. 
Alifanov and Mikhailov (1978) used the conjugate gradient method to solve the non-linear 
inverse thermal conductivity problem. Chen and Lin (1981) developed a pulse-spectrum tech- 
nique (PST) based iterative numerical algorithm for remote sensing of the thermal conduc- 
tivity of a non-homogcncous material for the one-dimensional (ID) case.lt is found that PST 
does give excellent results and is more robust in solving the inverse problem of the diffusion 
equation than that of wave equation.Jarny et al. (1986) employed an Output Least Square 
Method (OLSM) to determine thermal conductivities of materials. 

Tervola (1989) discussed the use of finite element method in conjunction with Davidson- 
Fletcher-Powell method for determination of thermal conductivity of a homogeneous mate- 
rial from the measured temperatures at certain ( finitely many ) points.The work concerns the 
situation where the thermal conductivity is dependent on temperature. Jarny et al. (1991) em- 
ployed a general optimization method using adjoint equation for solving multi-dimensional 
inverse heat conduction problems. Chen and co-workers (1996) employed the hybrid scheme 
of the Laplace Transform technique and the central difference approximation to estimate 
temperature dependent thermal conductivity from temperature measurements inside mate- 
rial at an arbitrary time. Lin and Chang (1997) proposed a numerical algorithm to estimate 
the thermal conductivity of a homogeneous material from boundary temperature measure- 
ments. Dowding et al. (1996) presented a laboratory method to measure the thermal prop- 
erties of a carbon-carbon composite material that is characterized by an orthotropic thermal 
conductivity and isotropic volumetric heat capacity. Alencar Jr. et al. (1998) presented a gen- 
eral solution for two-dimensional boundary inverse heat conduction problem by using the 
conjugate gradient method of minimization together with an elliptic scheme of numerical 
grid generation. Simulated measurements are used to illustrate the application of the present 
approach to the solution of an inverse problem of practical interest. 

The literature on the inverse determination of thermal conductivity from measured 
temperature data by finite-difference method are only a few. Lam and Yeung (1995) em- 
ployed a first order finite-difference method to determine thermal conductivity in a one- 
dimensional (ID) heat conduction domain using Cartesian coordinates. In a subsequent 
paper, Yeung and Lam (1996) have successfully extended their earlier results, examining the 
feasibility of using a second-order finite-difference technique to determine the thermal con- 
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ductivity in a one-dimensional (ID) heat conduction domain. Recently, Ghoshdastidar and 
Ray (1998) presented a second-order finite-difference procedure for the inverse determina- 
tion of thermal conductivity in two-dimensional (2D) square and cylindrical domains using 
available transient temperature data at discrete grid points. The procedure is capable of pre- 
dicting constant, spatial and temperature-dependent thermal conductivities. For 2D square 
and cylindrical domains, the estimated thermal conductivity is verified for a constant con- 
ductivity material such as steel using numerically simulated transient temperature data.For 
2D square domains, the results have also been tested against five bench mark solutions for 
ID domain in which thermal conductivities are either constant or linear/non-linear function 
of space and temperature by forcing a 2D domain to behave like a ID domain.The close 
agreement between the current results and the exact solutions confirm accuracy and effec- 
tiveness of the proposed finite-difference technique. It may be noted that the advantage of 
using finite-difference method is that the conductivity function can be obtained by solving 
the system of linear equations arising out of the discretization of the governing partial dif- 
ferential equation. Therefore, no prior information is required on the functional form of the 
thermal conductivity. 


1.3 Objectives of the present study 

The objectives of the present study are as follows. 

1. To extend the work of Ghoshdastidar and Ray (1998) to the special case of a 2D square 
domain of a known constant conductivity material (mild steel) in which no heat is gen- 
erated.Thus, the transient temperature data are created by obtaining numerically sim- 
ulated temperature for the situation when the body is heated to a uniform temperature 
(Tj) and then is allowed to cool by losing heat to the ambient. 

2. To develop an experimental set-up which will produce transient temperature data that 
will be used as the input to obtain the thermal conductivity of mild steel in a 2D square 
domain. 

3. Since, inverse heat conduction problems are very sensitive to measurement errors, an 
effort is made to quantify a realistic error in conductivity by adding random errors to 
the simulated temperature inputs ( 'T cxac t )■ 

4. To calculate conductivity using filtered T actua i where T ac t U ai = T eIQCt 4-measurement 
error 

5. To do a stability analysis based on 3 and 4. 
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1.4 Organization of the Thesis 


The present thesis has been organized in the following manner. 

> In Chapter 2, the finite-difference technique of Ghoshdastidar and Ray (1998) is dis- 
cussed with respect to the present problem. 

> Chapter 3, details the experimental set-up used to obtain transient temperature data 
required for the numerical algorithm and gives estimated conductivity values of mild 
steel using actual and filtered temperatures. The sources of experimental errors and 
possible ways to alleviate them are also mentioned. 

> Chapter 4 presents the results of the stability analysis based on numerically simulated 
temperature inputs ( T aclua i ) with the addition of the random measurement errors as 
well as that based on T actua i after filtering using Gram Orthogonal polynomial method 
described in Al-Khalidy (1998). 

> Finally in Chapter 5, the conclusion drawn from the study and the scope for the future 
work are indicated. 
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Chapter 2 


Problem Formulation 


2.1 Problem Statement . 

Consider the problem 'of unsteady heat conduction in a square plate (L x L) which 
is initially at a temperature, T t . Suddenly, all four sides are exposed to a cool ambient 
at Too. The problem can be assumed to be a two-dimensional one as the plate is thin so 
that the thermal picture in the plane normal to the Z-direction is more or less identical. The 
physical problem and computational domain is pictorially depicted in Fig.2. 1 . The simulated 
temperature profiles at various times have been computed using Explicit Finite-Difference 
scheme ( Appendix ’A’ ). 


2.2 Governing Differential Equation 


The governing differential equation for this case is, 

n r 21 — h ( d 2 T | d 2 T\ ■ dkdT ■ dkdT 
dt ~ K \ dx 2 ^ 8y 2 J " t " dxdx ' i ~ Oy dy 


( 2 . 1 ) 


Giving the spacewise as well as time wise temperature gradients as the input to the eq. 
(2.1) will result in the conductivity value ’k’ as output, ’k’ may or may not be uniform. The 
next article shows the details of the finite-difference scheme. 


2.3 Discretization of GDE 

Equation (2.1) is discretized at each grid point (i,j) and at present time, ’p’. In this 
case, all the space derivatives have been discretized using suitable second order difference 
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schemes. A square grid of Ax = Ay has been used. See Fig 2.2 


2.3.1 Lower-Left Corner Point (1,1) 
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2.3.5 Upper-Left Corner Point (1, n) 
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2.3.6 Interior Grid Points : i=2,n-l and j=2, n-1 
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2.3.7 Right Boundary : i=n and j=2, n-1 
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2.3.8 Top Boundary : i=2, n-1 and j=n 
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2.3.9 Upper-Right Corner Point (n, n) 
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2.4 Method of Solution 


Gauss-Seidel iterative method with successive under-relaxation was used to solve the 
set of linear simultaneous algebraic equations arising out of the discretization of eq.(2.1) at 
each grid point in the computational domain. The solution gives the conductivity value, ’k’ 
at each grid point. The accuracy of the solution depends on grid spacings, time increment 
and the times at which the temperature profiles are taken. The values of p and C v used 
for obtaining simulated temperatures are 7801 kg/m 3 and 473 J/kg.K respectively.lt was 
found that a minimum of 16 grid points are required to get a reasonably accurate solution 
for this case since the material ( Mild steel ) used in this study has the conductivity of 45 
W/m.K which was a prior known.The double-precision arithmetic was used. The rate of 
convergence was quite fast. The G-S iterative procedure was programmed in Fortran 77. The 
computations were performed using HP-9000 mainframe computer system of I.I.T. Kanpur. 
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Chapter 3 


Experimental Techniques and Results 


Experimental facilities and techniques used for the transient temperature measure- 
ments in a 4"x4"x0.125" specimen of mild steel to estimate its thermal conductivity are 
first described in this chapter.This is followed by the discussion of the results based on ac- 
tually measured data and after filtering the data.The reasons for the deviation from literature 
value of thermal conductivity of mild steel are given at the end of this chapter. 


3.1 Thermocouples 

A thermocouple is a small electric device used for measuring temperatures. It consists 
of two wires of dissimilar metals joined at their ends, forming two junctions. If the junctions 
are at different temperatures, a voltage will be produced that is approximately proportional 
to the temperature difference, a phenomenon discovered by the German physicist Thomas 
Seebeck. 

All temperature measurements in this study were taken with 26 SWG, Type J, insu- 
lated thermocouples, red one indicating iron and white one indicating constantan ( Fig.3.1. 
). According to Rangan et al. (1997), Type J thermocouples are made of iron as the positive 
thermocouple element versus constantan ( 57 % Cr, 43% Ni ) as the negative thermocou- 
ple element. They are generally recommended for use over the temperature range of -40°C 
to 800°C. The accuracy is usually ±2% in the range of -40 to 375°C, which the same is 
±0.5%in the range of 375 to 800°C. The Type J thermocouples are inexpensive, mechan- 
ically strong but rapidly deteriorates above 6O0°C'.In the present experiment temperature 
never exceeded 200°C. 
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3.2 Formation of Hot junction 


The thermocouple hot junctions were made by flaming in the Glass blowing Workshop, 
I.I.T. Kanpur. Prior to that the thermocouple ends were cleaned with the help of emery 
paper.The quality of the junctions were determined by applying certain tension. 


3.3 Calibration of Thermocouples 

Temperature-measuring devices should be calibrated periodically to evaluate their per- 
formance. Such calibration can be carried out either by comparison with a standard device 
whose accuracy is known. This method has the advantage that the test can be performed at 
any desired temperature range of interest, in a bath. The attainable accuracy depends on the 
uniformity of the temperature distribution in the test bath and the quality of the reference 
standards used. The standard liquid-in-glass thermometer can be used as a reference device. 
For tests in the range of -100 to 600°C, it is preferable to use an electrically-heated liquid 
bath. Since silicone oil baths are suitable for use as bath liquid in the temperature range 
of 50 to 250°C, all the 16 thermocouples used in this experiment were calibrated using an 
electrical heated silicone oil bath with a mercury-in-glass thermometer as the reference de- 
vice. The calibration was done in the Refrigeration and Air-conditioning laboratory of I.I.T. 
Kanpur. 


3.4 Heater 

The square specimen of Mild Steel was heated by electric 220/230 Volts AC “Heating 
Mantles” having energy regulator in the range of 0-100 W. The glass fibres enclosed in 
aluminum housing provides high thermal efficiency on the inner surface. The mantles are 
fitted with ON/OFF switches. The maker of the heater is Electro-Mech Instruments, Madras. 
See Fig. 3.1 for the photograph of the heater. 


3.5 Preparation of the Sample and Assembly 

A mild steel specimen of 4"x4"was cut from a Luge sheet of 0.125" thickness in the 
Central Workshop, I.I.T.Kanpur. The top and bottom faces were smoothed by polishing and 
the sides were smoothed by grinding. Next, the 16 equally spaced grid points were marked 
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on the top face. At the designated points micro holes were made so that the thermocouple 
hot junctions can be attached using a suitable joining technique. The distance between each 
thermocouple is |inches. In the present case, after experimenting with various joining tech- 
niques such as use of high temperature cement, spot welding etc., it was found that brazing 
gave the best joint. The brazing was carried out in the TA202 Metallurgy laboratory of I.I.T. 
Kanpur. 

After the thermocouples were brazed to the designated grid points, the strength of the 
joints were tested by using certain tension, the next step was to connect the cold junctions 
to the connectors by screwing which in turn are attached to two selector switches ( Fig. 3. 2. 
) each capable of handling 10 thermocouples. Finally, the positive and negative ends of the 
last switch were connected to the corresponding ends of the digital temperature recorder. It 
was ensured that the thermocouple number and the number on the switches tally. 

The temperature recorder ( Fig. 3. 2 ) has the range of ambient to 200°C. Needless to 
say that the recorder has ambient compensation built into it. The maker of the recorder is 
BlueBell, Kanpur. It has an' LCD display with one digit after the decimal.The accuracy of 
the recorder according to the supplier manual is ±0.5%. This was also checked against a 
mercury-in-glass thermometer reading of ambient as well as boiling water temperature. 


3.6 General Experimental Set-up and Procedure 

The Mild steel square sample supported by four tiny props was placed on an iron screen 
which in turn was positioned on the heater ( Fig.3.3 ). The heating mantle was then turned 
on at 50W setting for 28 minutes, 70W setting for 15 minutes and 90W setting for 2 minutes 
till a maximum temperature of 186.0 °C was shown on the digital display of the temperature 
recorder. The overall heating time was 45 minutes. The heating was essentially by radiation. 
Table 3.1 at “0” time shows the temperature at the designated grid points at this point in time. 
Clearly, uniform temperature could not be achieved. This is in contrast with the assumption 
in the numerical simulation that the specimen is heated to a uniform temperature. 

The next stage is the cooling of the sample. To do this the hot body was held carefully 
by a pair of tongs and transferred to a cooler surrounding ( at a temperature of ll.2°C ). The 
sample was supported by four aluminum props ( with sharp upper ends ) as shown in Fig.3.4. 
During cooling, the temperature at each grid point was recorded every five minutes till 45 
minutes have elapsed. The temperature history during cooling is shown in Table 3.1. The 
entire experiment was carried out in the Heat Transfer Laboratory. Fig.3.5 and Fig.3.6 show 
the actual thermocouple numbering and corresponding grid point numbering using double 
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subscript notation. 


3.7 Results based on Actually Measured Temperatures 

Table 3.2 shows the conductivity matrices calculated based on the actually measured 
temperatures at successive times starting from time zero. The results clearly are unsatisfac- 
tory as negative conductivity values or not a number (NaN) are not acceptable. However, 
this is not unexpected as it is well-known that inverse problems are highly sensitive to mea- 
surement errors in the input data. Even, 5% measurement errors can produce 50% errors in 
the output. 


3.8 Results based on Filtered Temperatures 

As mentioned in Art.3.7., the error into the results is usually greater than the error into 
the input data and the solution may be oscillating. Therefore, the solution may be useless 
when the real data are used. In this work a Gram Orthogonal Polynomial method with a 
moving averaging filter window is used for smoothing the noisy data. The method is based 
on a least square approximation. For more details, see Appendix ’B’. Table 3.3 shows the 
transient temperatures after 1 1 point filtering. Table 3.4 shows the conductivity matrices 
based on the filtered data. Even after filtering, there is virtually no improvement in the results 
pointing to a very high sensitivity of the present numerical algorithm to measurement errors. 
The possible resons why even after filtering the results did not improve may be attributed to 
the fact that only timewise filtering was done. Since the present problem is 2D, spacewise 
filtering also should have been performed to get better results. 


3.9 Sources of Experimental Errors 

Since the present algorithm is very sensitive to measurement errors, it is essential that 
the temperature measurements are as accurate as possible. There are several reasons for not 
obtaining the exact temperatures. Some of them are listed below. 

1. In the theoretical formulation, two-dimensional heat conduction is assumed. In actual 
practice, perfect two-dimensionality is very difficult to achieve. 
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2. There are 16 thermocouples used to measure temperatures at various grid points. These 
sensors are quite long and also thin. So, it produces fin effect thus decreasing the 
accuracy of temperature measurement. 

3. The props on which the plate was heated and cooled were made of materials other than 
mild steel and therefore, actual temperatures were not obtained. 

4. The thermocouples were brazed to the mild steel plate. The filler material for brazing 
is a copper-alloy ( copper, zinc and tin ). The filler material may also have distorted 
the temperature field. 

5. The 16 thermocouple reading should have been taken at the same instant. In the present 
case, there was a time lag since it is not humanly possible to record all the temperatures 
at exactly the same instant. 


3.10 Possible ways to avoid Experimental errors 

To improve the experimental set-up, the following is suggested A PC-based transient 
temperature measurement should be taken with the aid of data acquisition card. Use of non- 
invasive technique such as thermal imaging should be looked into. Better heating and cooling 
arrangements are desirable. Instead of Type J, Type K ( Chromel-Alumel ) thermocouple can 
be used for higher accuracy. High precision temperature recorder should be used. 
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Fig. 3. 3 The Photograph showing the heating arrangement 


Fig. 3. 4 The Photograph showing the cooling arrangement 
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Table. 3.1 Actually Measured Transient Temperatures during ct 





Table. 3.2 Conductivity Matrices based on Actually Measured Temperatures 
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NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

1.02 ' 

-2.30 

54.74 

-71.94 

-14.77 

-17.07 

14.35 

30.79 

-32.29 

-120.14 

-321.57 

529.78 

-8.82 

-15.34 

-90.12 

146.20 

14.22 
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33.17 

34.91 

-36.78 
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13.86 
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-13.19 

-15.78 
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14.59 

12.49 

- 76.16 

39.31 

-28.77 

-93.51 

-.27 

53.42 

1.63 

-9.94 

50.56 

- 27.33 

15.28 

2.30 
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-822.64 
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NaN 

NaN 
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-8.67 

-20.01 

7.73 

1.68 

-11.48 

-13.64 

8.25 

11.18 

11.26 

7.26 

-12.85 

-16.66 

4.99 

V 

13.31 

-5.54 

8.03 

J 




Table.3 3 Transient Temperatures during cooling after filtering 
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Table. 3 . ^Conductivity Matrices based on Filtered Temperatures 
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Chapter 4 


The Effect of Noise: Simulation Studies 


4.1 Introduction 

The previous chapter has dealt with determination of thermal conductivity of a ho- 
mogeneous plate using actual temperature measurements at discrete points and an inverse 
thermal conduction algorithm. It is evident that the estimated values of thermal conductivity 
are quite sensitive to the accuracy of the input temperature values. The errors at differ- 
ent points and at different times are, in general, not correlated. This can possibly lead to 
sharp spatial and temporal gradients in the temperature field, thereby affecting the calculated 
values of thermal conductivity. This results in severe inconsistencies such as prediction of 
negative or inaccurate conductivity values and spatially varying (in a random manner) con- 
ductivity values even when the test specimen is homogeneous. It is therefore essential to 
understand and quantify the effect of noise on such predictions to make such methods useful 
from a practical point of view. However, in an experimental set-up, it is difficult to control 
the level of noise as an independent parameter. An alternative approach is to simulate such 
problems as outlined in Figure 4.1. Initially, transient temperature data may be generated for 
a given problem (such as a plate of a specified material subjected to a certain temperature 
and then allowed to cool in atmosphere) using the direct heat-conduction formulation. Noise 
may be added to such data using randomly generated numbers within specified ranges. The 
resulting data may then be analysed using inverse calculation methods to predict the thermal 
conductivity of the material. The results can subsequently be compared with the conduc- 
tivity value initially used for the direct problem. If the inverse method being used is based 
on discretization of space, such as a finite difference technique, we generally obtain a range 
of k values for the discrete nodal points. In this case we need to compare the goodness of 
the solution based on some objective measure (parameter) which considers the range of k 
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values obtained. In this work, we propose a number of such criteria and evaluate their rela- 
tive merit. A simulation study also allows us to ascertain the utility of various filters before 
applying inverse methods to noisy data. In the following sections, we present the details of 
the simulation procedure employed and the inferences drawn. 


4.2 Problem Definition: Input Parameters 

We consider a homogeneous plate (4"x 4" x 0.125" ) made of steel. The plate is ini- 
tially at temperature 500°C and is suddenly exposed to ambient at time t=0. The properties 
of mild steel are specified as: 


Property 

Units 

Value 

Thermal conductivity, k 

W/mK 

45.0 

Density, p 

Kg/m 3 

7801.0 

Specific Heat, C p 

kJ/kgK 

473.0 


We calculate theoretical transient temperature data at 16 grid points. The points are chosen 
corresponding to the locations of thermocouple attachment as described in Chapter 3. 


4.3 Direct Heat-Conduction Formulation 

The direct heat-conduction formulation is employed to determine the transient temper- 
ature data. We assume that there is no temperature variation in the plate-thickness direction 
and accordingly solve the 2D heat-conduction equation: 

dT _ d 2 T \ 

dt a \ dx 2 dy 2 ) 

The boundary conditions and the method of solution using a finite difference formulation are 
described in Appendix ’A’. The step-size in time-increment is choosen in accordance with 
stability requirements ( of the numerical algorithm adopted ) described in Appendix A . 


4.4 Simulation of Noise 

In inverse heat conduction problem there are a number of measured quantities such 
as temperature, time, senser location, and specimen thickness. We assume that the various 
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quantities are accurately known except the temperature. The temperature measurements are 
assumed to contain the major sources of error or uncertainty. Any known systematic effects 
due to calibration errors, presence of sensors, conduction and convection losses or others 
are assumed to be removed to the extent that the remaining errors may be considered to be 
random. These random errors can then be statistically described. 

The direct problem, as described in Section 4.3 yields simulated temperature data for 
the inverse problem. Noise can then be added to this data in terms of random numbers within 
specified ranges. Thus, actual temperature T at a grid point may be given by 

T = T m + T' 

where T rn is the mean (exact) temperature and T' is the noise (measurement error). If maxi- 
mum noise level is specified as <5 , we obtain 

T — T m + r.S 

where 7' is a random number between -1 and +1. 

The algorithms available for generating random numbers generally yield pseudo ran- 
dom numbers. It is therefore important to consider the statistical description of the numbers 
generated by any such routine. 

One thousand such numbers were generated using NAG routine in the range -1 to +1. 
Their distribution is shown in Fig.4.2 and mean value and Standard deviation are 0.009 & 
0.582 respectively. 


4.5 Results and discussion 

It is of interest to consider the effect of noise in simulated temperatures on inverse cal- 
culation of k value. Since we are using a finite difference algorithm, the solution is obtained 
in terms of spatial distribution of k values.The goodness of solution may then be quantified 
on the basis of prediction at a given grid-point, or in a more average sense, in terms of a 
parameter which takes into account the entire field solution. We adopt the latter approach 
and consider various alternative parameters. 


4.5.1 Stability Parameters 

An initial input to the direct problem is the thermal conductivity of the plate material. 
For the 4x4 grid being used in the present case, 
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kth kth k t h k t h 

7 _ kth kth kth kth 

^ input — 

kth kth kth 

kth kth kth kth 

where k t h is the actual (theoretical) thermal conductivity of the plate material. The corre- 
sponding distribution of k values obtained on solving the inverse problem is 


kn ku kiz ku 

_ &21 &22 &23 &24 

K output — 

&31 ^32 &33 &34 

&41 ki2 ^43 &44 

where % (i=l-4. j=l-4) represents k values obtained at various grid point. The k output 
values will range between k in i n and k tnax . The following parameters are now suggested for 
quantifying the goodness of the results. 


Zl = 


kmin “f - k. 


max 


2% 


Z2 = 


kmax 

z . 

"'min 


Z 3 = 

24 : 24.1 = -^, 24.2 = ^» 

kth kth 

In the ideal case i.c.k max = k min = k th , Zl, Z2 and Z4 are equal to one and Z3=0. Zl & Z2 
should be considered together to avoid misleading interpretation under special circumstances 
such as k m(iX — k^/i T /\k and A'min — - ktti £^k Or kmin — ~Z~ ktti> In practice Zl alone 

without considering Z2 gave reasonable measure of the goodness of solution.i.e. such special 
case was not encountered in any of the simulations. 


( kmax kth) ( kmin kth ) 

kth 
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The parameter Z3 is a more objective measure of the goodness of the solution since 
it considers any deviation in least-square sense. Still Z1,Z2 and Z4 are considered in this 
work to demonstrate that physically unacceptable (negative ) values may be obtained if noise 
levels are high. The use of parameter Z3 fails to bring out this point. 

4.5.2 Truncation and Round-off errors 

The use of a numerical algorithm ( both direct and inverse problems ) involves trunca- 
tion and rounding-off errors. Therefore, in general k th and k output will not be identical even 
if no noise were added to the simulated temperatures. A series of typical calculations and the 
resulting Z values for time steps nl(nl=5, tl=l 16.68 sec) and n2 (n2=35, t2=816.76 sec) are 
shown in Table 4.1. It is evident that the values of these parameters deviate from the ideal 
values because of truncation and rounding-off errors.The deviation is typically 3% 

4.5.3 The effect of noise 

Now we investigate the effect of noise on the stability parameters. Similar to Table 4.1, 
we go through a series of calculations except that now we add random noise to the simulated 
temperatures. Since the temporal gradients change with time, we repeat this simulation for 
various times ( ranging from early to late times ). 

We are interested in the following : 

1 . What is the maximum noise level which still gives us reasonable estimate of k value 
(say, within ±10%of k t h) based on inverse calculations. 

2. At what stage of cooling is the numerical scheme most tolerant to noise, i.e., whether 
early_time data is preferred over late_time data. The temporal gradients of temperature 
are much steeper during early time compared to late_times. 

In order to study the maximum permissable noise level, the simulation was run a num- 
ber of times with S as a parameter. The value of S ranged between 0 and 1 (insteps of 0.01). 
Each simulation requires temperature information at two time instants (time steps nl and 
n2), which was generated using the direct formulation. Simulations were run for various 
pairs of nl and n2 values to encompass a range of temporal-gradient values including early 
and late times. 

Z values were then calculated for each simulation based on the k ou t pu t values obtained. 
Typical results on the effect of noise on k ou t pu t values and Z parameters as a function of 
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early/late time and maximum noise level (actual noise at a grid point = random number in 
the range -1 to +1 x maximum noise) are shown in Table 4,2.1 - 4.2.2, Fig. 4.3. In these 
figures, early time refers to n 1=5 and n2 =35 (i.e. tl =116.68 sec & t2 =816.76 sec) and late 
time refers to nl =875 and n2 =905 (i.e. tl =20418.97 sec k t2 = 21 1 19.05 sec).Note that 
the first row of temperature data in these tables corresponds to 5 = 0 , i.e. no noise has been 
added. The following trends emerge from comparision of these results. 

1 . From Fig. 4. 3 it is clear that the solution is close to theoretical k values up to a critical 
value of 6, and beyond it the solution is unstable often predicting physically unac- 
ceptable (negative) values. Z parameters can therefore be considered as parameters 
measuring the stability of solution. This critical value can be measured in terms of 
acceptable deviation in Z value. So, for example, we consider the inverse prediction 
to be acceptable if k is within 10% of k th . This gives Zl (acep(atlle ) as 0.9 - 1.1 and 
Mere rptabic) as 0.1414. Fig. 4.4 shows typically the determination of critical value 
of 5 based on Z3 parameter. This value will in general, be different at different times 
during the thermal history. Table 4.3 provides the critical 6 values ( S max ) based on Z1 
parameter for different times. 

2. Inverse calculations based on early-time data yield better results than late-time data; 
i.e., the numerical scheme can tolerate larger errors at early times. This is presumably 
related to steeper temporal gradients at early times. 

3. On calculating (based on simulations) the maximum permissable error in temperature 
at various times during the cooling of the specimen, we find ( Fig. 4.5. (b) ) a systematic 
variation of S max with t. This variation is described well in terms of an exponential 
decay function (Table 4.4). In (2) above, it is supposed that early time data yields better 
results presumably due to steeper gradients. To test this hypothesis, it was decided to 
plot S mnx w.r.t. The calculation of 5 max with The Values of ^ at various times 
were calculated using noise-free data. Fig. 4. 6 shows this variation, based both on Z1 
and Z3 parameters. It is found that a linear function (Table 4.5) describes the above 
data quite well, and that 6 rnnx has larger value for larger ^r. 

4. The simulations were also run for hypothetical materials with different values of k ( 
k =30, 45, 60, 75 ; same p and C p ). In each case, the results were described well 
by exponential decay function (S maz versus t, Fig. 4. 5) and linear function (S mQX versus 
f> Fig-4.7) 


That such a correlation is obtained based on numerical simulation studies, merits fur- 
ther attention;it is of interest to explain these results in terms of a theoretical frame work. 


27 



In ihe following section, we attempt to explain these trends in terms of a lumped parameter 
analysis. 


4.5.4 Lumped parameter analysis 

Consider a lumped body, P, suddenly exposed to a surrounding media, S (Fig.4.8), at 
a temperature T^. The body P loses the heat to the media S by conduction and the media S 
acts as a heat sink. We assume that the boundary between P and S is such that the equivalent 
heat transfer coefficient in the same as thermal conductivity of P. 

The governing equation for change in temperature of body P is 


pCpV 


dr 

dt 


kA 


(Too - T) 
L 


(4.1) 


where L is a characteristic length scale and is of the order j where V is the volume and A is 
the surface area of body P. For an initial temperature T 0 and with both k and independent 
of time, the solution of the above equation is 


T = Too ~ (Voo - T 0 ) exp 


( H ) 

l pC p L*} 


(4.2) 


where T 0 is the temperature of body P at time t=0. 

In the presence of disturbance, T = T m + T' and T 0 = T om + T' 0 where subscript m 
refers to exact (mean) temperature and subscript (') represents the disturbance in temperature. 
Note that the term T ' represents the same physical quantity as 6. A different symbol is used 
here to distinguish between the notations for numerical simulation and theoretical (lumped 
parameter) analysis. 

The above equation used in conjunction with noisy data will give reasonable estimates 
of k only up to a critical noise level T' inax \ above this noise level, the numerical algorithm 
becomes unstable. Therefore T max may be identified as the maximum noise level which is 
acceptable for describing the solution in terms of equation 4.2. 


Therefore we have, 


T a-T' 

- 1 m i -*■ max 


= Too - 


(Tqo (T 0 rn H - T 0 _ max )) exp 



(4.3) 
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also, for exact values of temperatures (i.e. if there were no noise), 

T,n = T„ - (Too - T om ) exp -4L.) . (4.4) 

subtracting eq.4.4 from eq.4.3 we obtain 


T — T 

max x r 


0 -max^P 



(4.5) 


The above equation implies that maximum permissable error in temperature measurement at 
a given time t which still gives reasonable ( acceptable ) values of k, decays exponentially 
with time.This is in agreement with Fig.4.5. Similar trend is obtained for different k th values 
(Fig.4.5). Also, from eq.4.5, 


T^;= exp {-^) 

Therefore a master plot between and kt can be obtained which should be independent 

1 o—max , 

of k^ value. Fig.4.9 shows such plot on normalised axes based on the results presented in 
Fig.4.5. 

As far as Fig. 4.10 is concerned, we can depict the theory behind it in the following 
manner. 

From equation (4.1), 


where a = jafe and 5 = L 

From equation ( 4.2 ), 


a 


dT 

dt 


= T-T O0 


(4.6) 


T - Too = ~ (Too - T 0 ) exp {-(dt) (4.7) 


where (3 = p- 
From eq. 4.6 & eq. 4.7 
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From equation ( 4.5 ) also 


(IT 

dt 


(4.8) 


~ ~~L? ^°° ~ To) ex P {-fit) 


T ' = T' 0 exp (-fit) (4.9) 

Dividing equation 4.9 by equation 4.8 , we obtain 


or 


t' t' t 2 

max ± p - Ly 

§ “ (T» -T„)'~a 


( Tp-mas ) ( l±PCp ) dT 

v k A^o-Too);' dt 


(4.10) 


The above equation implies that T max should vary linearly with ^ (since the other 

terms in the above equation are constant and T' 0 _ max is the maximum permissible noise at 

t=0). This is in good agreement with Fig.4.6-4.7. Again, the data for different k values 

t' 

is expected to merge if the results plotted on a normalised x-axis representing - °y° ■ ■ 
Although the values tend to come closer, (Fig.4.10), some variation is present presumably 
because all the details of a 2D problem are not captured in a lumped parameter analysis. 



In this section we consider the effect of filtering on the random errors that we imposed 
on the simulated temperatures. Here we are using the Gram Orthogonal Polynomial method 
for smoothing the noisy data, which we are getting after the addition of the random errors to 
the temperatures.The Gram Orthogonal Polynomial method has been described in Appendix 
B.The main idea is to findout that whether we are getting any improved thermal conductiv- 
ity values at the same level of noise after applying the filtering tcchnique.For that wc have 
chosen two time instants i.e for early time and for late time. At that time instants we got 
the simulated temperatures ( direct calculation ), and then added the error randomly. We then 
applied the filtering technique to improve the noisy data. The details are tabulated in the 
Table4.6. Based on these results it appears that this filtering technique is not so effective 
with the random errors. 
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Temperatures 



Fig-4.1 outline of the problem 

















max imn 



EaflyTime 
Late Time * 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

Maximum Noise in temperature, 8 ( C C) 

(a) 



80 


60 


40 

S .a 

F J 

20 

II 

0 

a 

-20 


-40 


-60 


-80 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Maximum Noise in temperature, 8 ( C C) 


Fig. 4.3 The Variation of parameters Zl,Z2,Z3and Z4 with maximum noise 
in temperature (8 ) at early ( nl= 5, n2= 35 ) and late ( nl=875, n2= 905 ) times 
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Fig.4.7 Gradient Vs 5 max (based on Z1 values) for different 'k' values 
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Table.4. 1 The effect of truncation and round-off errors in inverse problems 


theoritcal 


p = 7801 kg/m 3 C p = 473 kJ/kgK 

Direct Simulation 

At t = 116.68 sec At t = 816.76 sec 

456.81 458.% 458.% 

458.% 461.19 461.19 

458.96 461.19 461.19 

456.81 458.96 458.96 


490.90 

493.22 

493.22 

490.90 

493.22 

495.61 

495.61 

493.22 

493.22 

495.61 

495.61 

493.22 

490.90 

493.22 

493.22 

490.90 


45 

45 

45 

45 

45 

45 

45 

45 

45 

45 

45 

45 

45 

45 

45 

45 


Z1 = 1 
Z2 = 1 
Z3 = 0 

Z4: Z41 = 1 ; Z42 = 


No noise is added 


Inverse Calculation 







Q VO C4 
^ O ~T 


'A ^ O 
S p oc 
* ro O 



— r 4 r**, -r 

M M N M 


ro -r t 
N N N 


— • C J rO -r -»r 

N N N N N 


c-* m ro r-~ 

QO NO vO OO 
ro rb rb rb 
~T T NT T 


n O' On ro 

\C cq rq vq 
ro ro rb rb 
rr -T rr rr 

cn O' On ro 
p rq rq vO 

rb rb rb rb 
rr rr rr rr 


r*'* ro on r > -« 

OC vO no OO 
ro rb rb rb 
rr rr xr rr 


« s ? 


vo rt ; 1"' 

rb rb O 


ro — • 
On nO nO I""~ 


vo vO wr> — • 
co vo vo O 


rr <N OO 
vO NO CO 


ro 

CN 

vo 

■*£ 

ON 

OO 

O' 

, , 

* ‘ 

vo 

oo 

00 

p 

CN 

CN 

rq 

oi 

ZT 

oo 

NO 

oo 

vo 

00 

d 

*T 

rf 

ro 

ro 

r- 

VO 

p- 

rr 

NO 

'T 

— 

ro 

r< 

VO 


•o - 

p 

OO 


CN 

vq 

O 

00 

3 

ro 

ON 

— 


vb 

00 

r~"‘ 


vo 

r r 

rr 

ro 

rr 

rr 

wo 

vo 

OO 

r- 

ro 

ro 

rf 


or 

sO 

00 

o 

p ; 

ON 



3 

00 

CN 

o 

NO 

VO 

d 

Q 

rr 

ro 

VO 

vo 

rr 

rr 

vo 


r0 

ro 

OO 

rr 


rr 

ro 


On 

CN 

ON 

oo 

ON 

ro 

00 

CO 

r- 

CN 

oo 

rr 

o 

VO 

NO* 

vo 

P 

4 

WO* 

ro 

a 

ro 

ro 


rr 

8 

8 

OO 

On 

On 

vb 

OO 

OO 

wo 

VO 

WO 

rr 

rr 

rr 

8 

CN 

On 

On 

— * 

p 

oo 


—i 

wo 

vO 

NO 

*T 

-T 

rr 

nO 

vC 

O' 

O' 

On 

• — ■ 

— 

ob 

— 

— 

vo 

v O 

NO 

rT 

rr 

rr 

rr 

8 

8 

oc 

ON 

On 

nO 

DC 

oo 


rr rr rr rr 


VO VO VO vo 
rr rr rr rr 


rr rr TT T 


r- 

o 

3 

NO 

vo 

r- 

00 

vo 

oo 

p 

r- 

8 

On O 
p O 

CN 

s 

NO 

o 

ON* 

oo 

r-* 

r- 

OO ON 

NO* 

wo 

VO 

vo 

vo 

Wi 

vo 

VO VO 

VO 

rr 

rr 

rr 

rr 

rr 

rr 

rr rr 

rr 

ON 

NO 

vo 

r- 

ON 

ro 

CN NO 

s 

vo 

r- 

• 

'O 

Ol 

cc 

— 

p 

p 

p 

CN 

p 

p 

oo oo 

p 

On’ 

oo* 

— 


ON 

oo* 

d d 

On' 

wo 

vo 

NO 

NO 

vo 

VO 

$ 5? 

VO 

' r T 

r r 

rr 

rr 

rr 

rr 

rf 

OO 

rf 

vo 


Cl 

rr 

ro nQ 

_ 



Cl 

rr 

O 

ro 

O nO 

VO 

On 

p 

On 

ON 

p 

p 

rq vq 

vq 

00 

CO* 

o’ 

o 

OO 

On 

— — 

OO 

vo 

vo 

NO 

NO 

vo 

VO 

NO NO 

VO 

rf 

rr 

rr 

rr 

rr 

rr 

rr rr 

rr 

p» 

ro 

fN 

04 

Cl 


NO T 

rf 

~t 


o 

NO 

r- 

CN 

ro O 

r- 

p 

— 


oo 

P 

— 

p p 

vo 

NO 

d 

On 

oo 

r~ 

r- 

00 On 

NO 

VO 

vo 

VO 

wo 

vo 

vo 

VO vo 

VO 

"T 

r r 

rr 

rr 

rr 

r r 

rr T 

rr 


r~- 

p- 

oo 



O' 

n 

n 

OO 

rb 

O' 

ro 

O' 

S 

■^r 

rr 

rr 

CN 

CN 


NO 

NO 

CN 

vo 

vo 

rb 

ON 

ON 

On 

rr 

*T 

r r 

CN 

CN 

r- 

vq 

vq 

°i 

vo 

vb 

rb 

O' 

O' 

O' 

rr 

rr 

rr 

r^ 

r- 

oo 

•— « 

•— > • 

On 

C4 


OO 

rb 

O' 

ro 

O' 

8 

*r 

"T 

rr 


<— • O ro ON 
O' rO vo ~t 

o — ; r| 30 

— * rb cb d 
O' O' O' O' 
t rr rr rr 


On vn CN 

r- rr vo 
— • vo vo — « 
cb vo vo rO 
O' On On On 


vo VO CO *— 
co vq vo O 
rb vo vo ro 

O' O' O' On 

rr ^ rr rr 


rr OO O' O' 
CN O 'O *0 
O' ' r T O — * 


oo r- oo 
rs ro oo co 
CN rr VO vo 
ro vo vo ro 
O' O' O' O' 


VO r" — ro 
o rr NO vo 
— * ro <o On 
ro vo vo ri 
O' O' O' O' 
"T *t T ''T 


ro ro nO 

On VO — VO 

— ro ro — * 
O' O' O' 


— 04 OO \0 

2 3 a s 


rr rf r- ro 
P- fO On vo 
ON CN CN NO 

ri vo vo ro 

O' ON O' O' 


vo rr OO ro 
OO CN OO O 

ro p^ on oo 
o o o ri 

O' O' ON O' 


3 00 o OO 
00 P VO 
CN — Tf VO 


♦V’ 





Tabic 4.2.1 The effect of noise level on Z parameters for calrly lime steps 
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Table 4.2.2 The effect of noise level on Z parameters for late time steps 
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Table 4.3 Maximum permissible values & gradient for different times based on Z1 parameter 


S.No. 

Time 

nl 

Steps 

n2 

Gradient 

range 

Maximum 

Permissible 

Error 

1 

5 

35 

0.048685 

_ 

0.049168 

0.38- 

2 

35 

65 

0.045183 

- 

0.045638 

0.335 

3 

65 

95 

0.041938 

- 

0.04236 

0.275 

4 

95 

125 

0.038925 

- 

0.039317 

0.265 

5 

125 

155 

0.036129 

- 

0.036493 

0.29 

6 

155 

185 

0.033534 

- 

0.033872 

0.27 

7 

185 

215 

0.031125 

- 

0.031439 

0.245 

8 

215 

245 

0.02889 

- 

0.029181 

0.195 

9 

245 

275 

0.026814 

- 

0.027085 

0.15 

10 

275 

305 

0.024889 

- 

0.025139 

0.185 

11 

305 

335 

0.023101 

- 

0.023334 

0.17 

12 

335 

365 

0.021442 

- 

0.021658 

0.135 

13 

365 

395 

0.019902 

- 

0.020102 

0.13 

14 

395 

425 

0.018472 

- 

0.018658 

0.125 

15 

425 

455 

0.017145 

- 

0.017318 

0.105 

16 

455 

485 

0.015914 

- 

0.016074 

0.13 

17 

485 

515 

0.014771 

- 

0.014919 

0.1 

18 

515 

545 

0.01371 

- 

0.013848 

0.095 

19 

545 

575 

0.012725 

- 

0.012853 

0.07 

20 

575 

605 

0.011811 

- 

0.01193 

0.065 

21 

605 

635 

0.010962 

- 

0.011073 

0.085 

22 

635 

665 

0.010175 

- 

0.010277 

0.068 

23 

665 

695 

0.009444 

- 

0.009539 

0.058 

24 

695 

725 

0.008766 

- 

0.008854 

0.052 

25 

725 

755 

0.008136 


0.008218 

0.042 

26 

755 

785 

0.007552 

- 

0.007628 

0.047 

27 

785 

815 

0.007009 

- 

0.00708 

0.036 

28 

815 

845 

0.006506 

- 

0.006571 

0.034 

29 

845 

875 

0.006038 

- 

0.006099 

0.026 

30 

875 

905 

0.005605 

— 

0.005661 

0.024 


42 




Table 4.4 : The various terms describing the exponential equation 


Equation : 

y = y +Acxp(-Bt) 




k 

A 

B 

y o 

A 

30 

0.52647 

1.1 17e-04 

0 

0.0002619 

45 

0.39492 

1.106e-04 

-0.01174 

0.0002699 

60 

0.26526 

9.355e-05 

-0.01618 

3.747e-05 

75 

0.13425 

2.147e-04 

0.05592 

0.0001322 

Table 4.5 

Equation : 

: The various terms describing the best fit linear equation 

y = A+Bx 

k 

A 

B 

Mean 

S.D. 

30 

-0.0124 

10.4701 

0.99459 

0.01418 

45 

-0.0172 

7.73875 

0.98733 

0.01567 

60 

1.0I5c-5 

4.91502 

0.99462 

0.006 

75 

-0.00532 

3.61778 

0.94081 

0.01099 
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Table 4.6 Effect of Filtering on the Erroneous temperatures 





Chapter 5 


Conclusions and Scope for Future work 


This work presents a finite-difference based numerical algorithm for the inverse deter- 
mination of thermal conductivity in a 2D square domain using actually measured transient 
temperature data. It is found that for mild steel ( a constant conductivity material ) the algo- 
rithm is highly sensitive to the measurement errors in the input data, a fact not unexpected 
in inverse problems which are basically ill-posed. In addition to the experimental work, 
a detailed sensitivity analysis is also performed in order to find the maximum permissible 
measurement error ( in temperature ) for such problems. Interestingly, the study reveals that 
the maximum permissible measurement error is greater at early times and decreases expo- 
nentially with time. This means that early time temperature profiles should be used as input 
data. However, even at early times maximum allowable error in the temperature measure- 
ment is quite small which explains why the present experimental input data did not give 
correct conductivity value for mild steel. 

Future effort should be directed towards improving the experimental set-up so that 
high precision temperature measurement is possible. A PC-based experimentation through 
data acquisition card is highly recommended. A non-invasive temperature measurement 
technique such as thermal imaging can also be looked into. A better heating and cooling 
arrangement is necessary. 
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APPENDIX A 


Determination of Simulated Temperature profile 

using Finite-Difference 

A.l Governing Differential Equation 

The governing equation for the constant conductivty case is given by 

PC P f = * (0 + 0) (A.l) 

A.2 Discretization of GDE 

The details of the discretization of eq.(A.l) at each grid point is described next. See 
Fig.2.1 for the physical problem and computational domain and Fig. 2.2 for grid point num- 
bering. The difference schemes of first and second order accuracies are given in Appendix 
’C. 


A.2.1 Botom Boundary : j=l and i=2, n-1 


This surface is convecting heat to the surroundings.Therefore, 

= h(T - Too) 

h(rp rp } 

^ 2A y ~ k ^ 1 ' J 

=s- Tij-i = Tij+i - 2A y’i (T,j - T„) 
Discretizing equation (A.l), we get 


c nr-n., = k [ | Sn: 

A t K | (Ax) 2 ^ 


(Ay)‘ 


(A.2) 


A.2.2 Top Boundary : i=2, n-1 and j=n 

Tliis surface is convecting heat to the surroundings. Therefore, 
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■*f = h(T-T„) 


I i J 4 1 — 7 , 


2Ay — * (-^j 

=P T iJ+ 1 = T-^, - 2Ay| (Tij - r„) 

Discretizing equation (A.l), we get 

r ^i_ r .. _ r rP j _ 1 -2rf J .+j? i _ 1 -2Ay^(iy J -r 00 ) ) 

^ G P At { (Ax) z + (Ay) 2 J 


(A.3) 


A.2.3 Left Boundary : i=l and j=2, n-1 

This surface is convecting heat to the surroundings. Therefore, 

if = h(T - Too) 

h. (rp . rp ^ 

=* fc ~ fc l i ‘J 

=*• Ti-u- = r i+1 j - 2 Az£ (Tjj - Too) 

Discretizing equation (A.l), we get 

At = j + (AS? / 


(A.4) 


A.2.4 Right Boundary : i=l and j=2, n-1 

This surface is convecting heat to the surroundings. Therefore, 

~k % = h(T-T„) 

=>■ Ti+ij = 2i-ij “ ^Axjj; (Tij - ^oo) 
Discretizing equation (A.l), we get 


pC p ^ 


_ i j TUi- 2T rj+ T Ur^ [ T u-^ 4- > 

P At — K I (Ax) 2 J 


(A.5) 
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A.2.5 Interior Grid Points : i=2, n-1 and j=2, n-1 

Discretizing equation (A.l), we get 


rr-ip+1 

pC P T ^~ 


rpp 


At 


ii±L£ 


TLm 


(A i)' 


+ 


rpp 

i id ±l : 


(Avf 


A.2.6 Lower-Left Corner Point : (1,1) 


(A. 6) 


Discretizing equation (A.l), we get 


^ At <7^ + (arf 


(A. 7) 


A .2.7 Upper-Left Corner Point : (1, n) 


Discretizing equation (A.l), we get 


^ r ^ 1 -^_/J T H-i.>-2^+rf +lii -2Axi(rA-T tX) ) , I? J _ 1 -27f J +rf J _ 1 -2AyKrf J -r oe ) 
#,Gp At -*•( (S) 1 + <^7 — — 


(A.8) 


A.2.8 Lower-Right Corner Point : (n. 1) 


Discretizing equation (A.l), we get 


oC SSASi-t. I V-,.,-2V. i +tu,-2 

P^p At -fc | ^ 


(A.9) 


A.2.9 Upper-Right Corner Point : (n, n) 


Discretizing equation (A.l), we get 


rr nt'-n, 

P U PAi (A*) 3 + (A^ 


(A. 10) 


50 



CHECKING FOR STABILITY 

The coefficient of 7^ is taken from each of the equations ( A.2 ) - ( A. 10 ) and is seen that 
it must be greater or equal tozero. In the process, we get diffemet Ai's .To ensure stability 
we choose the minimum of the Ai's arising out of the stability condition. 
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Appendix B 

Digital Smoothing Filter : Gram Orthogonal 

Polynomial Method 

In this work a Gram orthogonal polynomial method ( Al-Khalidy, 1998 ) with a moving 
averaging filter windowsis used for smoothing the noisy data. This method is based on a least 
square approximation. This method does not need any information about the beginning and 
the end of the process. This is suitable for use in on-line methods of analysis. Application 
of digital filtering to a series of equally spaced seven data points is shown in fig.B.l. .Digital 
filter replaces each data value 1\ by a combination of itself and a number of adjacent nodes. 
Thus, 


fn(tk) = E bjPj(k) 

3 = 0 

where f(tk) is the smoothed value of the measured temperature at t=k. Let L is the number 
of points used to the left ( past temperatures ) and to the right ( future temperatures ) of the 
central point at time k. The subscript ’n’ refers to the total number of data points used for 
the smoothing process ( n=5 for L=2, n=7 for L=3, etc.). The parameters in eq. (B.l) can be 
calculated as 


, £,nk) Pj {k) 

°J ~ (2L+?+l)!(2L-j)! 
(2j+l)[(2L)\) 2 


Pj( k ) = 


J 

E 


m= 0 


(_l)™+j( ra + j)M(L + fc)N 
(m!) 2 (2£)M 


where T denotes the measured values of temperature. Thus, f(t) for each time step as the av- 
erage from T (t _i) to T (t+L) is calculated. This is sometimes called moving window average. 
Notice that = x(x — l)(x — 2) (x — m + 1) and m = 1, 2, 3, 

At the first steps the past temperatures can be set to be equal to the initial temperature. 
Similarly, at the last steps the future temperatures can be set equal to the final temperature. 
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The use of 7 to 1 1 data points is enough to obtain a good approximation. The temperature 
variation between two measurements should be large enough to see variation of temperature 
and should be greater than measurement errors. 


53 







Appendix C 

Difference Schemes Of First and Second Order 

Accuracy 


C.l Forward difference with error 0(Ax) 


7/ — Vi+l-Vi 

y i Ax 

— yi+2—2yi+l+yj 
(Ax ) 2 


C.2 Backward difference with error 0( Ax) 


7/ — Vi-Vi-l 

Vi Ax 

1," — Viz. 2 y»-i+y»- 2 

(Ax ) 2 


C.3 Forward difference with error 0(Ax) 2 


7/ — -y»+ 2 +%»+i- 3 yi 

y i 2Ax 

— y£+3+4yj+2— 5y,-+i +2y a - 

_ (Axp 


C.4 Backward difference with error O(Ax) 2 


' _ %-4yj_i+yj_2 
y i ~ 2Ax 

" _ 2y t — 5y,--i+4y,-_ 2 -y,-3 
* (Ax ) 2 


C.5 Central difference with error 0(Aa;) 2 


,/ _ yj+i-yj-i 
y i 2 Ax 


(c.l) 

(c.2) 


(c.3) 

(c.4) 


(c.5) 

(c.6) 


{c.l) 

(c.8) 


(c.9) 
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